Process and device for real-time spectral analysis of complex unsteady signals

ABSTRACT

The invention includes a real time spectral analysis of unsteady signals with complex values representatives of a physical phenomenon. Successive sections of a sample signal are modelled by a high order self-regressive process by estimating parameters characterizing the model by rapid Kalman filtering using an algorithm.

The invention relates to the real-time adaptive spectral analysis ofunstable signals having complex values representing a physicalphenomenon.

The invention can be used whenever it is necessary to analyse in realtime such signals which can be represented by a self-regressingparametric model, by assuming that the coefficients of the model remainsteady over short time intervals. It applies more particularly to allthe fields of medical imagery, for example using NMR, and more generallywhenever spectral analysis must be made on organic signals, for examplein electrocardiography and electro encephalography. It finds however aparticularly important application in the spectral analysis of thesignals furnished by ultrasonic Doppler velocimeters for studying theblood flow in the heart cavities or the vessels.

At the present time, parametric spectral analysis is relatively littleused in Doppler velocimetry and in medical ultrasonic echography. Itremains at the laboratory stage. It consists in modelling the signalsusing a parametric model, generally of the adjusted mean self-regressingtype, whose parameters are adjusted by optimization using algorithms ofthe least square type. This approach has the advantage of working onshort signal sections but it involves using low order models and hasconsequently limited frequency resolution, which is particularlytroublesome for the spectral analysis of very unsteady signals.Moreover, the known methods do not process the demodulated signals withcomplex values as such, but process each of the components separately,respectively in phase and in quadrature.

Commercial apparatuses use non parametric spectral analysis methods ofFourier analysis type so as to be able to use them with rapid Fouriertransformation algorithms. But these methods require working on longsections of the signals so as to obtain a suitable estimate of thedesired spectral representation, which is incompatible with the natureof the unsteady states met with in the above mentioned fields.

The invention aims more particularly at providing a spectral analysisprocess and device answering better than those known heretofore therequirements of practice, particularly in that they improve thefrequency resolution, make possible a more accurate analysis of thesignals and restrict the variability of the spectra obtained forrepresenting flows in closely related conditions.

To this end, the invention proposes more particularly an adaptivespectral analysis process by modelling the phenomenon represented by thesignal characterized in that successive sections of the sampled signalare modelled, using a high order self-regressing process, close or equalto the maximum authorized by the duration of the sections and thesampling frequency, by making an estimate of the parametercharacterizing the model by rapid Kalman filtering using an algorithmwhich is simplified in the full limit authorized by the local steadystate of the model in the section and in that the spectral power densityis determined from the estimated parameters; to avoid risks ofinstability related to the choice of the high order, the estimate of theparameters of the self-regressing model are regularized by using aconstraint which consists in minimizing a criterion of the form

    J.sub.2 (a)=J.sup.1 (a)+Ω(a)

where J₁ (a) is the sum of the least squares [y(n)-z(n)]², z designatingthe sampled function modelling the signal via the parameters a to bedetermined and y the function represented by the samples of the signaland where Ω(a) is a regularizing function taking into account the apriori knowledge of the nature of the function or obtained byoptimization.

Through using a high order model, very good frequency definition isobtained. The process makes it possible to obtain parameters directlyfor calculating the spectral power density of the signal in the analysiswindow. It is in particular possible to use a model order equal to thenumber of sampling points in the Doppler signal, which gives anexcellent frequency definition, whereas the present-day spectralanalysis methods using a self-regressing model and criterion of theleast squares have an order of size three to seven times less than thenumber of sampled points.

The invention also provides a spectral analysis device which can bedirectly fitted to a conventional Doppler velocimeter, without modifyingthe existing probe and its associated electronic circuits, nor itsdisplay means.

The invention will be better understood from the following descriptionof the process and device according to the invention and the reasoningwhich led to it. It refers to the accompanying drawings in which:

FIG. 1 is a general diagram showing the parameters used in Dopplervelocimetry of the blood flow in the particular case of a continuousDoppler technique;

FIG. 2 is a diagram showing the envelope of a representative Dopplerspectrum;

FIG. 3 is a general diagram showing demodulation of the Doppler signal;and

FIGS. 4 and 5 are block diagrams showing a possible implementation ofthe process of the invention.

MODELLING THE UNSTEADY SIGNALS

Before starting the description of the invention, information will berecalled concerning the modelling of the signals representing unsteadyphenomena and concerning the regularization of the model.

In the case where the unsteady state of the signal is sufficiently slowwith respect to the speed of variation of the amplitude of the signal,the latter may be considered as "locally steady", i.e. steady over ananalysis section of N samples, N being typically a few tens or so. Thesamples observed y(n), n=1,2,. . .,N are described by a self-regressingmodel (AR) expressed by:

    y(n)=z(n)+b(n)                                             (1) ##EQU1## where a.sup.t is the vector of the coefficients of the model whose order is p and where b(n) represents the procedure generating the signal, of unknown variance σ.sup.2.

The spectral analysis is then made following a criterion which is oftenthat of the least squares, which consists in minimizing J₁ (a): ##EQU2##

Since the model is linear, the calculating algorithms are rapid. Butobtaining a good approximation is subordinated to the choice of an orderp as high as possible; now, if we take p close to N or equal thereto,the minimization of J₁ is not a satisfactory criterion because of theexcessive variance of the solution.

To overcome this difficulty and stabilize the solution, it has beenproposed to regularize the model by adding a stabilizing function to J₁,which makes it possible to keep a high value of p (close to or evenequal to N). The criterion then consists in minimizing J₂ (a):

    J.sub.2 (a)=J.sub.1 (a)+Ω(a)                         (5)

where (a) is a regularizing function , numerous examples of which exist.

This function will however as a general rule be of the form:

    Ω(a)=|a-a.sub.o |.sup.t M(a-a.sub.o)

Ωis therefore scalar, product of the transposed matrix |a-a_(O) |^(t),of a positive symmetrical and defined matrix of dimension p.p (p beingthe dimension of a) and (a-a_(O)). The transpose is obviously a linevector.

If M is the identity matrix, we arrive at the simplified form which willbe given at (7) further on.

The choice of the regularizing function Ω (a) translates the a prioriinformation about the nature of the local spectral properties of theanalyzed signal. This information is that the unknown spectral powerdensity is a "smooth" function, i.e. having a certain regularity,compared with an unorganized spectrum formed of a multitude of spectralrays without any relation with each other. For reasons of convenience ofimplementation, a quadratic function is chosen , which has the advantageof leading to an explicit and linear solution with respect to the valuesobserved y(n).

A measure of "spectral smoothness" is formed by ##EQU3## where H(f) isthe frequency transfer function of the bleaching filter for the analyzedsignal, which is also the inverse of its generating filter which isdescribed by the equations (1) and (2). A high value of D_(k) signifiesa non smooth spectral power density, all the more so the higher theorder k of the drift. A zero order spectral smoothness is used: ##EQU4##which is, except for a constant, the euclidean norm of the vector of theparameters a, which leads in most cases to adopting the regularizingfunctional: ##EQU5## where the coefficient μ is a positive real numbercalled "regularizing coefficient" and which, as will be seen, plays avery great role in the operation of the invention and which is chosen byan a posteriori likelihood maximum method.

REGULARIZATION MODELLING ACCORDING TO THE INVENTION

The invention uses a procedure which may be considered as in two steps,based particularly on the fact that, in practice, the signals are veryoften signals for which a local study of the phenomenon may be made bydividing the time horizon into blocks of length N and that an assumptionmay be made of a local steady state in each block, so as to be able touse the properties of local invariance of the model by then carrying outan adaptive spectral analysis reflecting a slow variation of theparameters with respect to the amplitude fluctuations of the signal,namely:

1. A locally steady model is sought, i.e. an equivalent invariant modeldefined by the optimum of the criterion (5) calculated over a reducedtime window, by making a compromise between too great a window length(which does not allow the unsteady states of the signal to be followedfaithfully) and too short a length (which does not allow the model to beestimated under good conditions), which leads to choosing an AR model ofthe highest order p possible, close to N or equal to N, which involvesstabilizing the model.

This model is sought for the value of the regularization coefficient andthat of the variance σ² of the generator process which make thelikelihood of the model maximum.

2. The calculations for maximizing the likelihood of the model aresimplified by adopting a sub-optimum solution:

the variation range of μ of the equation (6) of the

minimization criterion is reduced to a finite preselected set ofdiscrete values, and

using an algorithm which may be rapid, the optimum solution a andpossibly the corresponding likelihood V(σ², μ/y) are calculated fromsamples of the observed signal, so as to select the highest.

The first step leads to organizing the samples of the signal inconsecutive blocks of length N and carrying out regularization on eachblock. A constraint is imposed relating the set of parameters locallyand translating spectral smoothness information: that leads tominimizing, on each of the blocks, a criterion of the form:

    J.sub.2 (a)=J.sub.1 (a)+μ|a-a.sub.O |.sup.2(7)

When the process generating the analyzed signal is assumed Gaussian,that is tantamount to admitting that the parameters a are characterized,in a Bayesian interpretation, by the law of a priori probability:##EQU6##

The vector of the AR parameters is initialized by choosing a_(O=) 0 inthe case of a steady signal. It will be seen further on that this mannerof proceeding may be kept for certain unsteady signals. In the othercases, for each current block, the result of calculating the precedingblock will be chosen when the blocks are adjacent.

It still remains to fix the value of μ and σ² for minimizing thecriterion (7). This choice is made, as was mentioned above, bymaximizing the likelihood of the model. It will be seen that μ and σ²,called hyperparameters, may be decoupled in this maximization.

To explain the first step, the equations of state associated with theproblem should be given and it should be shown how the standard Kalmanfilter (which is based on a Riccati equation not using the local steadystate of the model which is at the basis of the adaptive methods, northe property of shift of the successive vectors y_(p) (n) which is theconsequence of the choice of an AR model) a simpler filter may besubstituted, using Chandrasekhar's equations for calculating theparameters a which minimize, with μ and σ² fixed, the criterion (7), aswell as for calculating the corresponding likelihood. It will be seenthat this filter, with invariant state model, allows the regularizationconstraint to be introduced.

For that, these properties are caused to appear explicitly by definingan "extended vector-parameter" or state vector with m components:

    a.sub.m (i)=[0.sup.t i-1, .sup.a p,.sup.a p-1,. . ., a.sub.1 1,0,0,. . .,0].sup.t, with m>p

and the overall vector of the data is considered:

    y.sub.m =y(p+1),y(-p+2),. . ., y(0), y(1),. . .,y(m-p).sup.t.

The linear model defined by the equations (1) and (2) may then be in theform of the following equivalent state model:

    a.sub.m (i+1)=Da.sub.m (i)i=1, 2, . . .

    y(i)=y.sup.t.sub.m a.sub.m (i)+b(i)

where D is the shift operator: ##EQU7##

The calculation of the solution may then be effected by the followingKalman filter:

    a .sub.m (i+1/i)=Da.sub.m (i/i-1)+k.sub.m (i)r(i).sup.-1 y(i)-y.sup.t.sub.m a.sub.m (i/i-1)                                           (8) ##EQU8## with the initial conditions:  where y.sub.m designates the simply conjugated vector of vector y.sub.m and k.sup.* designates the conjugate transposed matrix of k.

The introduction of such additional variables is apparentlyunfavourable, for it increases the dimensions of the state vector. Butthis impression is misleading for such introduction causes theinvariance properties to appear explicitly, which are used.

It may moreover be verified that, at each recurrence, the number of nonzero coordinates of the gain vector k_(m) (i) remains equal to p.

Since this state model is invariant (D and y_(m) being constant) andwith steady variance (σ² being constant over a block), a rapid algorithmis obtained directly by replacing the standard equations byChandrasekhar's equations. Not only do these equations allow rapidresolution, but they further make it possible to pass naturally tosquare root algorithms offering good digital stability.

Hereafter, by "algorithm B-CAR" will be designated the regularizationrules using the algorithm, applied to blocks of a size appropriate tothe nature of the signal. These equations only propagate the incrementsof the nominal quantities of the equation (9), which are themselvesfactorized. The direct application of such factorization techniques tothe filter (8) (9) leads to the following equations: ##EQU9##

In these equations, the fact is used that M and P_(m) are matrices withHermitian symmetry: M^(*) =M and P^(*) _(m) =P_(m). Condensed writing ofthe algorithm is obtained by projection, by only causing the componentsto appear which effectively come into play at each recurrence in thecalculation. For that, an effective gain-vector k(i) is defined, ofdimension p, containing the non zero coordinates of k_(m) (i):

    k.sub.m (i)=[0.sup.t.sub.i+1, k(i).sup.t, 0.sup.t.sub.m-p-i-l ].sup.t

Similarly, if α is the dimension of M(i), we may also write:

    B(i)=[0.sub.α,k,S(i).sup.t, 0.sub.α,m-p-i-l ].sup.t

where S(i) is a matrix of dimensions [(p+1),α]. The algorithm may thenbe in the form: ##EQU10## the updating of the parameters taking place inaccordance with:

    a((i+1/i)+a(i/i-l)+k(i)r(i).sup.-1 [y()-y .sub.p (i).sup.t a(e/e-l)](13)

This algorithm has the advantage of applying not only to "pre-windowed"problems (processing of the steady case of an isolated block in theunsteady case) but also to covariance type problems (processing ofadjacent blocks in the unsteady case).

Initialization of the algorithm is conditioned by the a prioricovariance matrix P (1/0) from which B(1) and M(1) are derived, at thetime of factorization of the initial increment:

    δP.sub.m (2)=P.sub.M (2:1)-P.sub.m (1/0)=B(1)M(1)B(1) .sup.t

The initial covariance increment may also be written, returning to theequations, in extended form:

    δP.sub.m (2)=DP.sub.m (1/0) D.sup.t -k.sub.m (1)r(1).sup.-1 k.sub.m (1) .sup.* -P.sub.m (1/0

with ##STR1##

In the general case, if we assume ##EQU11## then we may resolve##EQU12## where ##EQU13##

The matrix P(1/0) is chosen diagonal and the rank of δ_(m) P(2) isgenerally δ=3. Since the initial factorization is not unique, thesignature matrix of P_(m) (2) is advantageously chosen for M(1), whichgives the beginning of a square root algorithm whose numerical qualitiesare better.

If the a priori covariance matrix P(1/0) is diagonal, the vectors gm anddm are identically zero and the initial increment P_(m) (2) is reducedto:

    δP.sub.m (2)=-v.sub.0 S.sub.m S.sup.*.sub.m +ν.sub.0 v.sub.m v.sup.*.sub.m -k.sub.m (1r(1).sup.-1 k.sub.m (1.sup.*,

of a rank at most equal to 3 since it is broken down into a sum of threedyadic or antiscalar matrices, each of a rank equal to 1. δP_(m) (2) mayalso be written in the form:

    δP.sub.m (2)=B(1)M(1)B(1).sup.*

with, as central factor= ##EQU14## which leads, once the projection hasbeen made so as obtain the reduced dimension algorithm, to: ##EQU15##under the conditions, which are generally fulfilled, that:

the coefficients a_(i) of the model, which are complex, have real andimaginary parts whose a priori laws of probability are normal, centred,independent and of variance (1/2)τ² (we thus find again the a prioricovariance matrix P(1/0)=τ² I),

the generator noise of the observed signal b(n), which is also complex,has real and imaginary parts whose a priori laws of probability are alsonormal, centred, independent and of variance (1/2)σ².

This rigorous solution uses the values effectively observed prior to theinitial time of the window in the general vector y_(m). This method,which may be termed "covariance", assumes operating:

either on an isolated data block of a length twice the chosen order p,

or on a data block adjacent a prior block whose contents have been kept.

The down-count of the elementary arithmetic operations, limited tomultiplications alone, shows that the complexity of this algorithm isθ(11p) by recurrence. If we choose p=N, it can be seen that the totalcomplexity of processing N data is θ(N²), which is normal for analgorithm of this type.

A solution for simplifying the algorithm may be termed "pre-windowing".It assumes that the p values prior to the initial time in themeasurement window are zero, with for consequences:

    k.sub.m (1)=D P.sub.m (1/0)y.sub.m =0 ##EQU16##

    δP.sub.m (2)=v.sub.0 s.sub.m s.sup.*.sub.m +v.sub.O v.sub.m v.sup.* m

The rank of the initial increment isα=2. As initial factor we maychoose: ##EQU17##

The complexity of the calculation is substantially reduced.

When the nature of the phenomenon represented by the signal is wellknown, the hyperparameters may be chosen a priori. In the opposite case,calculation of the likelihood of the hyperparameters is useful.

With a Kalman filter, the likelihood of the hyperparameters (or itslogarithm) may be calculated by recurrence, using quantities interveningin the algorithm.

The search for the likelihood maximum amounts in fact to a problem ofregularized least squares already represented by the equation (7). Thisproblem has a simple Bayesian interpretation since the minimization of(7) is equivalent to maximizing the likelihood criterion V(a): ##EQU18##

When the hyperparameters μ, σ² are used, the calculation shows that thesolution a is in fact the mean of the a posteriori law defined by theconditional distribution of the data: ##EQU19## and the a priori law:##EQU20##

We may then consider that this a priori law specifies a class ofestimators via the parameters μ and σ². These are the hyperparameters ofthe problem. Since the laws are normal, the likelihood of thehyperparameters is represented by an integral:

    V(σ.sup.2,μ,y)=∫f(y/a,σ.sup.2)f(a/μ,σ.sup.2)da(16)

The model to be chosen is that which maximizes this likelihood withrespect to the hyperparameters. In this sense, the method is adaptivesince the choice of a priori values may depend on the data.

Finally, by sequential application of Bayes rule to the relation, weobtain: ##EQU21##

The conditional marginal density of y(n) is used which is possible toobtain from the observed samples y(1), y(2),. . ., y(n-1): ##EQU22##

The calculation of this conditional density requires knowingσ² and a.But,

    y(n)=y.sub.p (n).sup.t a(n/n-1)+e(n)

where a(n/n-1) is the optimum estimate (i.e. conditional to the past(y(1), y(2),. . .,y(n-1)), and where e(n) is the innovation of theobserved process. We also then have: ##EQU23##

The variance of the innovationσ_(e) (n)² remains to be evaluated. But aKalman filter is only sensitive to the relative variance alone of theobservation noise (i.e. here of the generator process) and of the statenoise (i.e. here of the a priori variance of the parameters). Since therecurrent calculation of a(n/n-1) only therefore depends on μ and not onthe value of σ², we may take arbitrarily σ² =1 in the algorithm: we thenhave:

    σ.sub.e (n).sup.2 r(n)=σ.sup.2 [y.sup.t.sub.m P.sub.m (n/n-1)y.sub.m +1]

whence we derive: ##EQU24##

This formula shows that the two hyperparameters are decoupled. Theestimate σ² of σ² according to the likelihood maximum (marginal) isthen: ##EQU25## and the (marginal likelihood of μ is then: ##EQU26##

To sum up, the complete process comprises, when we are to determine theoptimum hyperparameters:

for each of the discrete values of μ chosen for the calculation of J₂,the estimation of a vector a(N/N) with the algorithm (10(11), with σ²=1,

the choice of the parameter a and of the hyperparameters σ² and μcorresponding to the highest likelihood; and

calculation of the corresponding spectral power density g(f) using theconventional relation: ##EQU27## for

    -1/2≦f≦+1/2

APPLICATION OF THE MODELLING ACCORDING TO THE INVENTION TO DOPPLER BLOODVELOCIMETRY SIGNALS

The choice of the duration of the self regressive model windows orblocks, of the number of samples to be taken per block, of the order ofthe model and of the regularizing criterion involves knowing theprincipal characters of the Doppler signals. These signals are obtainedby a circuit of the kind shown schematically in FIG. 1. The speed fieldstudied is illuminated by a monochromatic or pulsed acoustic wave comingfrom a transmitter 10 and the wave is measured which is reflected by thematerial targets (red blood corpuscles) moving under the action of thisspeed field, collected by a receiver 12 comprising a demodulator.

If the speed field were steady and uniform, the signal received r(t)could be derived from the transmitted signal a(t) by a simple frequencytranslation. This is not the case, because of the contribution of alarge number of elementary targets having a dispersion with respect tothe ideal mean target; the signal received has a continuous spectrum, ofa narrow band width with respect to the frequency of the carrier (a fewkHzas compared with a few MHz). In the frequent case of a pulsedtransmission Doppler signal, the form of the transmitted wave also comesinto play. In both cases a narrow band spectrum is obtained of the kindshown in FIG. 2, which is limited to two symmetrical bands of width 2Fcentred at -f₀ and f₀ respectively, F being very much less than thefrequency f₀ of the carrier.

To simplify the processing of the signals, existing apparatis usedemodulation, generally synchronous, which leads to breaking up thesignal received r(t) into two components, respectively in phase p(t) andin quadrature q(t). These low frequency (of a few kHz) signals require areduced sampling rate for their useful spectrum is reduced by a factorof about 1000 with respect to the carrier. The circuit may be the oneshown schematically in FIG. 3.

In a complex representation we may write: ##EQU28## where ω₀ is theangular frequency the carrier and where e(t) designates the complexmodulation envelope of the transmitted signal, and ##EQU29## where p(t)and q(t) are low frequency signals generated by the demodulator 14 ofFIG. 3, forming respectively the real part and the imaginary part of thecomplex value "Doppler signal" y(t).

Calculation shows that the spectrum of the signal received r(t) can beobtained by Fourier transformation of the covariance of the Dopplersignal y(t). Moreover, p(t) and q(t) are generally correlated, whichprohibits assuming--nevertheless admitted in usual processing of theDoppler signal--that the mutual covariance R_(PQ) (τ) of p and q ortheir interspectral power density G_(PQ) (f) is zero. Now, thisassumption is necessary if p(t) and q(t) are processed separately.

The invention therefore processes p(t) and q(t) simultaneously, which ispossible because the above described spectral analysis method has beenin fact designed to be applied to complex value signals.

CHOICE OF THE REGULATION PARAMETERS IN THE APPLICATION TO DOPPLERVELOCIMETRY

Three regulation parameters must be chosen, i.e. a priori (which will bethe general case since the nature of the Doppler signal is known), orfollowing an optimization procedure.

It is a question of:

N, duration of an analysis window or length of a data block,

p, order of the AR model adopted for describing the Doppler signal,

μ, regularization coefficient or hyperparameter.

The initialization method must also be chosen. It is defined by:

the choice of the initial parameters a₀ at the beginning of each block:a₀ may be either zero or else the result of processing of the precedingblock;

the choice between the so-called "covariance" version and that called"pre-windowed".

The choice of the three parameters N, p, and μ results from acompromise. To increase N and p improves the quality of the model, but alimitation is caused by the unsteadiness of the signal which leads toreducing N and by the calculation load which leads to reducing p.

Good results are obtained with a sampling frequency of some tens of kHzand with observation windows of a maximum length of N=64 points, i.e. amaximum duration of 1 to 2 ms (instead of 5 to 10 ms in present daycommercial apparatus). Since there is no need for a spectral analysis ofthe Doppler signal at such short time intervals, the sampled data may beorganized into blocks of N=64 points at most and only working on oneblock out of ten. Thus, a spectral power density can be obtained every10 to 20 ms at most, i.e. at the rate of about fifty per cardiac period.

The initialization may be achieved by pre-windowing or by covariance, asmentioned above. Pre-windowing, as used at the present time, reduces thenumber of multiplications to be made but adversely affects thelikelihood of the model, especially if p is close to N. The covariancemethod reduces the maximum order of the model to p=N/2 if therecurrences are begun at time p+1 (in the case of isolated windows) andit is a little more complex.

The choice of the hyperparameter μ is essential for it governs thedegree of smoothness of the solution and must maximize the likelihood ofthe model of the signal with respect to the data: the values of the ARparameters do not in fact depend on the variances σ² and τ² but solelyon their ratio μ. Experience shows that the curve or variation of thelikelihood as a function of μ varies very little close to the maximumand that the results of the analysis are little affected by variationsof an order of size of this parameter. Experience also shows that, for agiven Doppler signal, the variations of the optimum value of μ from onedata block to another are also of an order of size. The complexity ofthe method may therefore be reduced by substituting, in the step formaximizing the likelihood of the model which must be carried out on eachblock of data, a previous step for classifying the signals as a functionof the constant values to be attributed to

For this to be possible, it is indispensable to standardize μ withrespect to the power of the Doppler signal. In fact, the coefficientintervenes in the minimized criterion:

    J.sub.2 (a)=J.sub.1 (a+μ|a-a.sub.0 |.sup.2(7)

where J₁ (a) has the dimension of the energy of the signal and the ARparameters are without dimension. The coefficient μ has then a dimensionand varies as the square of the scale in which the Doppler signal ismeasured.

To make μ independent of the power of the signal, it is then necessaryto choose:

    μ=μσ.sup.2 (N/p)

where μ₀ is the new regulation parameter, this time without dimension,and σ² is the variance of the signal which may be furnished by theacquisition module of the device used (24). Such standardization makesit possible to work with a constant μ₀ throughout a cardiac cycle,independently of the power fluctuations of the Doppler signal.

The invention is susceptible of numerous embodiments. In particular, thesensitivity of the process can be reduced to the limited accuracy of thecalculations by propagating the square roots of the covariance matricesrather than the matrices themselves: such a modification is simple tocarry out, as mentioned above.

CONSTRUCTION OF THE DEVICE

The device used for implementing the process may have the generalconstruction shown in FIG. 4, which seeks a compromise between thequality of the results and the volume of calculations, is adaptable toexisting apparatus and uses the intrinsic possibilities of digitalsignal processors at present available.

The device comprises two processors coupled for calculating the Dopplerspectrum. The first one, 20, calculates the parameters of a model of theDoppler signal received using the algorithm B-CAR. The second one, 22,delivers the spectrum from this model. This second processor 22 controlsan acquisition card 24 which acquires the two Doppler channelssimultaneously, at a sampling frequency which is for example of 40 kHzfor studying the blood flow, chosen so as to optimize the compromisebetween the duration of the sections analyzed and the frequency accuracyof the spectral analysis. To avoid possible dissymmetry in theacquisition of the Doppler signals, a single converter is used. The cardmay be adapted for detecting possible saturation of the Doppler signals.

The processor 20 is programmed for scaling the hyperparameter μ of themethod from the estimate of the instantaneous power of the Dopplersignal delivered by processor 22, and calculating the AR parameters ofthe Doppler signal on each data block transmitted by the same module.

The processor 22 calculates the spectrum of the signal on a block, fromthe AR parameters delivered by processor 20 under the control of thedisplay module 26 which fixes the number of points on the axis of thefrequencies and the number of grey levels; it may be adapted forcalculating a possible DC component in the acquisition chain. It must inaddition evaluate the instantaneous power of the signal, used byprocessor 20 for optimizing the results delivered.

The architecture of the apparatus may be as shown in FIG. 5. Processor22 executes the algorithm for calculating the spectral power density andcontrols the data input and output interfaces 28-30-32. Processor 20executes the algorithm B-CAR.

A host mini computer forms a master machine which loads the programmeinto the processors and controls the execution of the programmes. A ROMimplantation of the programme would avoid the need for a host minicomputer. The acquisition 30 and display 32 interfaces transfer thedata, on the one hand between the acquisition device and processor 22and between the display device 26 and this same processor 22. Eachprocessor has its own working memory 36.

It is necessary to carry out the calculations in a 32 bit format withfloating point. Circuits available commercially respond to such needs.It is a question on the one hand of signal processors formed of a 32 bitcomputation unit with floating point, a control unit, memories andinterface devices integrated in a single circuit (for example DSP32 ofATT and μ/PD77230 of NEC), on the other hand, 32 bit arithmetic units.The process has in particular been put into practice using DSP32processors whose internal buses are accessible from outside the case(which makes it possible to increase the memory size by 4 kbytes to 56kbytes), and which have a parallel interface for DMA memory accesswhatever the activity of the processor and a high flow rate seriesinterface (8 Mbits/second).

We claim:
 1. Process for the real time spectral analysis of unsteadysignals with complex values, representing a physical phenomenon, inwhich the phenomenon represented by the signal is modelled by aself-regressive process, characterized in that:successive sections ofthe sampled signal are modelled by a high order self-regressive process,close to or equal to the maximum authorized by the duration of thesections and the sampling frequency, by estimating the parameterscharacterizing the model by rapid Kalman filtering using an algorithmwhich is simplified to the extent authorized by the local steady stateof the model in the section and in that the spectral power density g(f)is determined from the estimated parameters; and in that: the parametersa of the self-regressive model are regularized: ##EQU30## where a is thevector of the coefficients of the model whose order is p, where b(n)represents the process generating the signal, where z(n) designates thesamples of the function representing the phenomenon and where y(n)designates the observed samples, by using a constraint which consists inminimizing a criterion of the form:

    J.sub.2 (a)=J.sub.1 (a)+Ω(a)

where J₁ (a) is the sum of the least squares

    [y(n)-z(n)].sup.2,

and where Ω(a) is a regularizing function taking into account the apriori knowledge of the nature of the phenomenon.
 2. Process accordingto claim 1, characterized in that Ω(a) is of the form:

    Ω(a)=μ|a-a.sub.0 |.sup.2

where a₀ is a predetermined value, chosen a priori zero or formed by theresult obtained for the preceding section.
 3. Process according to claim2, characterized in that μ is taken equal to μσ² where μ is apredetermined value and σ² the estimated power of the signal.
 4. Processaccording to claim 2, characterized in that the spectral power densityis calculated by the relation: ##EQU31## .
 5. Process according to claim2, characterized in that the in phase and in quadrature components of ademodulated narrow band signal are processed simultaneously.
 6. Processaccording to claim 2, characterized in that, with the signals beingdemodulated Doppler signals representative of the blood flow having afrequency spectrum of a few kHz, sections of about 64 samples are usedtaken at a frequency of about 40 kHz and the in phase and in quadraturecomponents of the demodulated signal are processed simultaneously.